clear;
load('options_BE.mat')
load('Data_FR_emp_prof.mat', 'Data_FR_emp_prof')
nsect = length(Data_FR_emp_prof);
%nsect=1
kjInd_FR_emp_prof = zeros(nsect,1);

for s = 1:nsect
    sigma=Data_FR_emp_prof(s,3);    
    %sigma=1.0
    k_cool = @(k_j)sd_loge_int(k_j, sigma);
    x = fsolve(k_cool,5, options_BE)
    kjInd_FR_emp_prof(s,1)=x;
    save('kjInd_FR_emp_prof.mat', 'kjInd_FR_emp_prof');
end



save('kjInd_FR_emp_prof.mat', 'kjInd_FR_emp_prof');


kappa1 = zeros(nsect,1);
kappa2 = zeros(nsect,1);
theta = zeros(nsect,1);

for s = 1:nsect
fun3 = @(z, k) (1-(z.^2)).*(z.*(exp(z))).^k .*exp(z);
Integral3 = integral(@(z)fun3(z, kjInd_FR_emp_prof(s,1)),0,1);

kappa1(s,1)=kjInd_FR_emp_prof(s,1).*exp(-kjInd_FR_emp_prof(s,1)-1).*Integral3;
kappa2(s,1)=kappa1(s,1)./kjInd_FR_emp_prof(s,1);
theta(s,1)=(kappa2(s,1).*(kjInd_FR_emp_prof(s,1)+1))./((1./(kjInd_FR_emp_prof(s,1)+1))-(kappa1(s,1)+kappa2(s,1))  );
end
thetajInd_FR_emp_prof=theta;
save('thetajInd_FR_emp_prof.mat', 'thetajInd_FR_emp_prof');


init = ones(nsect,1);
CC=Data_FR_emp_prof(:,4);
AA=theta;
ns=s;
beta_cool= @(B)eqsystem(AA, B, CC, ns);
    x = fsolve(beta_cool, init, options_BE)
    betajInd_FR_emp_prof=x;
    save('betajInd_FR_emp_prof.mat', 'betajInd_FR_emp_prof');
    
    
intra_dist_FR_emp_prof = zeros(nsect,1);
intra_dist_FR_emp_prof=100.*( ((kappa2.*((kjInd_FR_emp_prof+1).^2)).^(-1./(kjInd_FR_emp_prof+1))) -1  );
trick=betajInd_FR_emp_prof'*thetajInd_FR_emp_prof;
inter_dist_FR_emp_prof = zeros(nsect,1);
inter_dist_FR_emp_prof=100.*((thetajInd_FR_emp_prof./(trick))-1);

        save('inter_dist_FR_emp_prof.mat', 'inter_dist_FR_emp_prof');
        save('intra_dist_FR_emp_prof.mat', 'intra_dist_FR_emp_prof');